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Abstract 


This  report  describes  a  numerical  method  that  may  be  used  to  efficiently  locate  and 
track  magnetic  targets  with  a  tensor  gradiometer.  A  target  containing  ferromagnetic  material 
can  be  adequately  modeled  at  a  distance  by  an  equivalent  magnetic  dipole.  This  magnetic 
target  can  be  observed  by  means  of  a  magnetic  gradiometer  that  measures  a  symmetric, 
traceless  gradient  tensor  as  a  function  of  time.  Of  interest  is  the  inverse  problem  of  the 
determination  of  the  magnetic  parameters  of  the  target,  and  its  position  and  velocity  relative  to 
the  sensor  at  each  time  step.  The  previous  method  of  direct  inversion  of  the  non-linear 
equations  of  the  magnetic  gradient  tensor  provided  multiple  solutions,  and  the  results  can  be 
highly  sensitive  to  noise  in  data. 

In  this  study,  the  determination  of  target  magnetic  moment,  position  and  velocity  is 
formulated  as  an  optimal  stochastic  estimation  problem,  which  could  be  solved  using  a 
sequential  Monte  Carlo  based  approach  known  as  the  ‘particle  filter’.  In  addition  to  the 
conventional  particle  filter,  the  proposed  tracking  and  classification  algorithm  uses  the 
unscented  Kalman  filter  (UKF)  to  generate  the  prior  distribution  of  the  unknown  parameters. 

The  proposed  method  is  then  demonstrated  by  using  it  to  locate  and  track  an 
automobile  over  a  period  of  time  using  real  data  collected  with  a  magnetic  gradiometer.  Two 
cases  are  investigated:  (i)  the  observation  contains  only  gradiometer  data  when  a  double 
solution  exists,  and  (ii)  magnetic  field  components  are  added  to  the  previous  case  and  a 
unique  solution  is  obtained.  The  automobile  was  moving  either  on  a  straight  or  a  curved  track. 


Resume 

Le  present  rapport  decrit  une  methode  numerique  qui  peut  etre  utilisee  pour  localiser  et 
poursuivre  avec  efficacite  des  cibles  magnetiques  au  moyen  d’un  gradiometre  de  tenseur.  Une 
cible  contenant  un  materiau  ferromagnetique  peut  etre  modelee  a  distance  de  fa9on 
satisfaisante  au  moyen  d’un  doublet  magnetique  equivalent.  Cette  cible  magnetique  peut  etre 
observee  a  l’aide  d’un  gradiometre  magnetique  qui  mesure  un  tenseur  de  gradient  symetrique 
en  fonction  du  temps.  Ce  qui  nous  interesse  c’est  le  probleme  inverse  de  la  determination  des 
parametres  magnetiques  de  la  cible  et  de  ses  position  et  vitesse  par  rapport  au  capteur  a 
chaque  intervalle  de  temps.  L’ancienne  methode  d’inversion  directe  des  equations  non 
lineaires  du  tenseur  de  gradient  magnetique  donnait  lieu  a  des  solutions  multiples,  et  les 
resultats  peuvent  etre  tres  sensibles  au  bruit  contenu  dans  les  donnees. 

Dans  la  presente  etude,  la  determination  du  moment  magnetique,  de  la  position  et  de  la 
vitesse  de  la  cible  est  formulee  comme  un  probleme  d’ estimation  stochastique  optimale,  qui 
pourrait  etre  resolu  a  l’aide  d’une  methode  sequentielle  de  Monte  Carlo  appelee  «  fibre 
particulaire ».  En  plus  du  filtre  particulaire  classique,  Talgorithme  de  poursuite  et  de 
classification  propose  fait  appel  a  Testimateur  de  Julier  et  Uhlmann  pour  generer  la 
distribution  a  priori  des  parametres  inconnus. 

La  methode  proposee  fait  ensuite  l’objet  d’une  demonstration  pour  localiser  et 
poursuivre  une  automobile  pendant  un  certain  temps  a  l’aide  de  donnees  reelles  collectees  au 
moyen  d’un  gradiometre  magnetique.  Deux  cas  sont  etudies  :  (i)  1’ observation  est  basee  sur 
des  donnees  de  gradiometre  seulement,  ce  qui  donne  lieu  a  deux  solutions;  (ii)  des 
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composantes  du  champ  magnetique  sont  ajoutees  au  cas  precedent,  ce  qui  produit  une  solution 
unique.  L’automobile  se  dcplagait  soit  sur  une  piste  rectiligne,  soit  sur  une  piste  courbe. 
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Executive  summary 


Introduction 

A  requirement  exists  to  improve  the  Navy’s  capability  to  conduct  intelligence, 
surveillance  and  reconnaissance  (ISR)  operations.  Precise  determination  of  target  motion 
parameters,  i.e.  position,  velocity,  and  target  classification,  are  primary  concerns  in  automated 
surveillance  systems.  A  moving  target  containing  ferromagnetic  material  can  be  observed  by 
means  of  a  magnetic  tensor  gradiometer  that  measure  a  symmetric  gradient  tensor  as  a 
function  of  time.  Of  interest  is  the  inverse  problem  of  the  determination  of  the  position  and 
magnetic  parameters  of  the  target  at  a  given  time  step  from  its  magnetic  signature  collected  up 
to  and  including  that  time  step.  The  previous  method  of  direct  inversion  of  the  non-linear 
equations  of  the  magnetic  gradient  tensor  provided  multiple  solutions  in  addition  to  the 
physical  one,  and  the  results  can  be  highly  sensitive  to  noise  in  the  data. 

Work  description  and  results 

This  report  describes  a  novel  numerical  method  that  may  be  used  to  efficiently  locate 
and  track  magnetic  targets  with  a  tensor  gradiometer.  The  determination  of  target  magnetic 
moment,  position  and  velocity  is  formulated  as  an  optimal  stochastic  estimation  problem, 
which  could  be  solved  using  a  sequential  Monte  Carlo  based  approach  known  as  the  ‘particle 
filter’.  In  addition  to  the  conventional  particle  filter,  the  proposed  tracking  and  classification 
algorithm  uses  the  unscented  Kalman  filter  (UKF)  to  generate  the  prior  distribution  of  the 
unknown  parameters. 

To  analyze  the  dynamic  system  of  a  moving  target  one  needs  two  models:  (i)  a  model 
describing  the  evolution  of  the  state  with  time,  and  (ii)  a  model  relating  the  noisy 
measurements  to  the  state.  It  is  demonstrated  that,  if  the  system  dynamics  is  imposed  (e.g. 
rectilinear  motion)  and  only  the  gradient  data  are  used,  one  obtains  the  physical  solution 
together  with  a  second  solution  representing  the  reflection  through  the  origin.  This  study  then 
shows  that  the  rotationally  invariant  quantities  associated  with  the  gradient  tensor  have 
powerful  properties  for  target  localization. 

Finally,  using  real  data  collected  with  a  magnetic  gradient  sensor,  the  ability  of  the 
algorithm  to  locate  and  track  an  automobile  moving  either  on  a  straight  or  a  curved  path  is 
demonstrated.  Two  cases  are  analyzed:  (i)  the  observation  contains  only  gradiometer  data 
when  a  double  solution  given  by  a  scaled  state  vector  exists,  and  (ii)  magnetic  field 
components  are  added  to  the  previous  case  and  all  target  parameters  are  uniquely  determined. 

Significance  and  future  work 

The  need  for  tracking  under  realistic  conditions  has  motivated  a  series  of  trials  where 
magnetic  signals  from  various  targets  were  recorded.  The  proposed  method  of  recursive 
Bayesian  estimation  proved  to  be  accurate  and  applicable  to  various  situations  (straight  or 
curved  paths).  The  validation  of  the  presented  algorithms  on  experimental  data  indicates  the 
possibility  to  incorporate  it  into  automated  surveillance  systems. 


Birsan  M.  2006.  Recursive  Bayesian  method  for  tracking  a  magnetic  target  with  a 
gradiometer.  DRDC  Atlantic  TM  2006-188.  Defense  R&D  Canada  -  Atlantic. 
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Sommaire 


Introduction 

II  existe  un  besoin  d’ameliorer  la  capacite  des  forces  navales  pour  mener  des 
operations  de  renseignement,  de  surveillance  et  de  reconnaissance  (RSR).  La  determination 
avec  precision  des  parametres  de  deplacement  de  la  cible,  c.-a-d.  la  position,  la  vitesse  et  la 
classification  de  la  cible,  est  une  preoccupation  principale  dans  les  systemes  de  surveillance 
automatises.  On  peut  observer  une  cible  mobile  contenant  du  materiau  ferromagnetique  au 
moyen  d’un  gradiometre  de  tenseur  magnetique  qui  mesure  un  tenseur  de  gradient  symetrique 
en  fonction  du  temps.  Ce  qui  nous  interesse  c’est  le  probleme  inverse  de  la  determination  de 
la  position  et  des  parametres  magnetiques  de  la  cible  a  un  intervalle  de  temps  donne  a  partir 
de  sa  signature  magnetique  collectee  jusqu’a  cet  intervalle  de  temps.  L’ancienne  methode 
d’ inversion  directe  des  equations  non  lineaires  du  tenseur  de  gradient  magnetique  donnait  lieu 
a  des  solutions  multiples  en  plus  de  la  solution  physique,  et  les  resultats  peuvent  etre  tres 
sensibles  au  bruit  contenu  dans  les  donnees. 

Description  des  recherches  et  resultats 

Le  present  rapport  decrit  une  methode  numerique  novatrice  qui  peut  etre  utilisee  pour 
localiser  et  poursuivre  avec  efficacite  des  cibles  magnetiques  au  moyen  d’un  gradiometre  de 
tenseur.  La  determination  du  moment  magnetique,  de  la  position  et  de  la  vitesse  de  la  cible  est 
formulee  comme  un  probleme  d’estimation  stochastique  optimale,  qui  pourrait  etre  resolu  a 
l’aide  d’une  methode  sequentielle  de  Monte  Carlo  appelee  «  filtre  particulaire  ».  En  plus  du 
filtre  particulaire  classique,  Lalgorithme  de  poursuite  et  de  classification  propose  fait  appel  a 
l’estimateur  de  Julier  et  Uhlmann  pour  generer  la  distribution  a  priori  des  parametres 
inconnus. 

Pour  analyser  le  systeme  dynamique  d’une  cible  mobile,  il  faut  deux  modeles  :  (i)  un 
modele  decrivant  revolution  de  l’etat  dans  le  temps;  (ii)  un  modele  mettant  en 
correspondance  les  mesures  bruitees  et  l’etat.  II  est  demontre  que,  si  la  dynamique  du  systeme 
est  imposee  (p.  ex.  un  mouvement  rectiligne)  et  que  seules  les  donnees  de  gradient  soient 
utilisees,  on  obtient  la  solution  physique  ainsi  qu’une  deuxieme  solution  representant  la 
reflexion  passant  par  l’origine.  Cette  etude  demontre  done  que  les  grandeurs  invariantes  en 
rotation  associees  au  tenseur  de  gradient  ont  des  proprietes  puissantes  pour  la  localisation  des 
cibles. 

Finalement,  on  fait  la  demonstration  de  la  capacite  de  Lalgorithme  a  localiser  et  a 
poursuivre,  a  l’aide  de  donnees  reelles  collectees  au  moyen  d’un  tenseur  de  gradient 
magnetique,  une  automobile  se  depla5ant  selon  une  trajectoire  rectiligne  ou  courbe.  Deux  cas 
sont  analyses  :  (i)  l’observation  est  basee  sur  des  donnees  de  gradiometre  seulement,  ce  qui 
donne  lieu  a  deux  solutions  obtenues  d’un  vecteur  d’etat  mis  a  l’echelle;  (ii)  des  composantes 
du  champ  magnetique  sont  ajoutees  au  cas  precedent,  et  tous  les  parametres  de  cible  sont 
determines  de  fag  on  univoque. 
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Portee  et  recherches  futures 


Le  besoin  d’effectuer  la  poursuite  dans  des  conditions  realistes  a  donne  lieu  a  une  serie 
d’essais  dans  lesquels  on  a  enregistre  les  signatures  magnetiques  de  diverses  cibles.  La 
methode  proposee  faisant  appel  a  la  solution  bayesienne  recursive  s’est  averee  exacte  et 
applicable  a  diverses  situations  (trajectoires  rectilignes  ou  courbes).  La  validation  des 
algorithmes  presentes  au  moyen  de  donnees  experimentales  indique  la  possibility  de  les 
incorporer  dans  des  systemes  de  surveillance  automatises. 
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1.  Introduction 


Precise  determination  of  target  motion  parameters,  i.e.  position,  velocity,  and  target 
classification,  are  primary  concerns  in  automated  surveillance  systems.  A  target  containing 
ferromagnetic  material  can  be  adequately  modeled  at  a  distance  by  an  equivalent  magnetic 
dipole  moment.  This  target  in  motion  can  be  observed  by  means  of  a  magnetic  tensor 
gradiometer  that  measure  a  symmetric,  traceless  gradient  tensor  as  a  function  of  time  as  it 
passes.  Of  interest  is  the  inverse  problem  of  the  determination  of  the  position  and  magnetic 
parameters  of  the  target  at  time  step  k  from  its  magnetic  signature  collected  up  to  and 
including  time  k. 

Wynn  [1]  demonstrated  by  numerical  methods  that,  using  the  measurements  from  one 
magnetic  gradiometer  only,  the  inversion  solution  is  given  by  two  scaled  vectors  representing 
the  source-to-sensor  bearing  vector,  and  the  magnetic  moment  orientation  vector.  The  range 
and  the  magnetic  moment  magnitude  cannot  be  resolved  with  this  method.  Moreover,  there 
were  4  solutions  for  bearing  vector  and  scaled  moment  vector,  with  one  being  the  physical 
solution,  the  second  being  called  the  ‘ghost’  solution,  and  two  more  solutions  obtained  by 
reflecting  the  first  two  through  the  origin  centered  on  the  measurement  point.  If  an  additional 
five  measured  quantities  are  provided  given  by  the  rate  of  change  of  the  gradient  tensor  [2,  8], 
the  solution  represented  by  the  three  scaled  vectors  (position,  velocity  and  moment)  is  unique 
except  for  the  reflection  through  the  origin.  However,  the  results  obtained  by  the  direct 
inversion  of  the  non-linear  equations  can  be  highly  sensitive  to  noise  in  data.  For  this  reason, 
novel  signal  processing  techniques  more  robust  against  noise  in  data  were  developed  based  on 
the  statistical  approach  of  finding  the  best  correlation  to  possible  signal.  Such  methods  are 
matched-field  processing,  linear  statistical  analysis,  and  Monte  Carlo  simulation. 

In  this  report,  we  will  concentrate  to  a  sequential  Monte  Carlo  method  for  target 
tracking  known  as  the  “particle  filter”.  The  approach  makes  use  of  the  recursive  Bayesian 
estimation  (filtering)  technique  and  of  the  assumption  that  the  states  follow  a  first  order 
Markov  process.  Thus,  difference  equations  are  used  to  model  the  evolution  of  the  system 
with  time  and  measurements  are  assumed  to  be  available  at  discrete  times.  Using  the  state- 
space  approach  to  modeling  dynamic  systems  and  the  discrete-time  formulation  of  the 
problem,  the  aim  is  to  estimate  the  hidden  state  process  from  the  measurements. 

Let  define  xk  as  the  state  of  the  system  at  time  step  k,  and  zi:k  =  {z\,  z2,  ...,  zk}  as  the 
observation  (measurements)  history  of  a  system  from  time  1  to  k.  The  state  variables  are  the 
position,  velocity  and  magnetic  moment  of  the  target.  Because  of  either  noise  in  the  state 
evolution  process  or  uncertainty  as  to  the  exact  nature  of  the  process  itself,  the  state  vector  xk 
is  generally  regarded  as  a  random  variable.  In  the  Bayesian  filtering  technique,  one  attempts 
to  construct  an  estimate  of  the  posterior  probability  density  function  (pdf),  p(xk  |  z1:k).  Since 
all  information  provided  by  zkk  is  conveyed  by  the  posterior  density,  it  may  be  said  to  be  the 
complete  solution  to  the  estimation  problem.  A  recursive  Bayesian  algorithm  imposes  the 
constrain  that  the  estimate  of  p(xk  |  zi:k)  should  be  generated  solely  from  the  previous  posterior 
density,  p(xk_i  |  z k _ ,  ) ,  and  the  most  recent  measurement  zk.  In  this  way,  it  is  not  necessary  to 
store  the  complete  data  set  or  to  reprocess  existing  data  when  a  new  measurement  becomes 
available. 

The  recursive  propagation  of  the  posterior  density  is  a  conceptual  solution  only  that 
can  be  determined  analytically  only  in  a  restrictive  set  of  cases.  When  the  analytical  solution 
is  intractable,  a  Monte  Carlo  based  approach  to  recursive  Bayesian  filtering,  called  the  particle 
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filter,  is  one  method  that  approximates  the  optimal  Bayesian  solution.  In  the  Monte  Carlo 
method,  a  set  of  random  samples  (particles)  are  drawn  from  a  target  distribution,  such  as 
p(x  |  z).  In  general,  this  distribution  is  not  known.  We  will  use  q(xk  |  z1:k)  *  p(xk  |  zi:k)  to 
denote  a  proposal  distribution  from  which  samples  can  be  drawn.  The  main  drawback  of  the 
conventional  particle  filter  is  that  it  uses  transition  prior,  p(xk  |  xk_i),  as  the  proposal 
distribution.  The  transition  prior  does  not  take  into  account  current  observation  data.  To 
overcome  this  difficulty,  the  unscented  Kalman  filter  (UKF)  was  proposed  to  generate  better 
proposal  distributions  by  taking  into  consideration  the  most  recent  observation. 

To  analyze  the  dynamic  system  of  a  moving  target,  one  needs  two  models:  (i)  a  model 
describing  the  evolution  of  the  state  with  time,  and  (ii)  a  model  relating  the  noisy 
measurements  to  the  state.  It  is  shown  in  this  study  that,  if  the  system  dynamics  is  imposed 
(e.g.  rectilinear  motion)  and  only  the  gradient  data  are  used,  one  obtains  the  same  solution  as 
the  one  obtained  by  Wynn  using  10  measurements  from  the  gradient  and  gradient-rate  tensors. 
If  the  magnetic  field  measurements  are  usable  simultaneously  with  the  gradient  data,  we 
obtain  three  additional  measurement  points  for  a  total  of  eight,  allowing  us  to  resolve  both  the 
range-moment  ambiguity  and  the  ghost  problem. 

In  this  report,  the  theory  of  target  localization  with  a  gradiometer-magnetometer  array 
is  first  reviewed  and  the  target  localization  properties  of  the  rotation  invariant  quantities 
associated  with  the  gradient  tensor  are  investigated.  Then,  using  real  data  collected  by  a 
gradiometer  sensor,  the  ability  of  the  recursive  Bayesian  algorithm  to  track  an  automobile 
over  a  period  of  time  is  demonstrated.  Two  cases  are  analyzed:  (i)  the  observation  contains 
only  gradiometer  data  when  a  double  solution  given  by  the  scaled  state  vector  exists,  and  (ii) 
magnetic  field  components  are  added  to  the  previous  case  when  all  target  parameters  are 
uniquely  determined.  It  will  be  shown  that  the  operational  range  of  the  gradiometer  is  limited 
by  the  actual  noise  level  that  depends  on  both  the  gradiometer  construction  and  the  system 
noise. 
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2.  Locating  a  magnetic  target  with  a  gradiometer 


The  mathematical  model  used  for  the  target  is  a  moving  magnetic  dipole.  The  target  is 
fully  characterized  by  its  position  and  the  value  of  the  magnetic  dipole  moment.  Let  consider 
that  r  is  the  position  vector  from  the  source  to  the  point  of  observation  in  meters,  m  is  the 
magnetic  moment  vector  of  the  dipole,  and  V  is  the  velocity  vector  in  m/sec.  The  magnetic 
flux  density  vector  B  at  a  given  point  due  to  a  magnetic  dipole  is  given  by  the  formula: 


4  n 


3  (r,m)r 


I  |5 

r 


(1) 


where  p(>  is  the  permeability  of  the  air  (=  4ti1  0"7),  and  <•,•>  is  the  inner  product  operator.  In 
the  above  formula,  all  vectors  are  defined  in  the  same  coordinate  system,  which  normally  is 
the  sensor  coordinates.  The  three  spatial  derivatives  of  the  three  components  of  the  magnetic 
flux  density  vector  B  generate  the  magnetic  field  gradient  tensor: 


SBj  _  Mo 
dr,  4  n 


15  (r,m  )rirj  3{mjrj  +  mjri  + 


(2) 


where  8y  is  the  Kronecker’s  delta.  There  are  only  five  independent  components  of  the 
gradient  tensor  (2)  because  V  •  B  and  V  x  B  are  equal  to  zero  in  source-free  lossless  regions. 
However,  this  is  substantially  more  than  the  three  data  points  provided  by  a  vector 
magnetometer. 

Equation  (2)  was  inverted  by  Wynn  [1]  to  obtain  five  parameters  for  the  scaled 
moment  vector,  M  =  3//0m/4/zv4 ,  and  source-to-sensor  bearing  (orientation)  vector, 
R  =  r/r ,  in  terms  of  the  components  of  the  gradient  tensor.  One  orientation  parameter  is  not 
independent  because  the  sum  of  the  squared  directional  cosines  equals  one.  Of  course,  the 
sixth  unknown  is  the  length  of  the  vector  r  that  cannot  be  resolved.  This  is  called  the  “range- 
moment  ambiguity”.  A  further  complication  arises  from  the  fact  that  the  equations  to  be 
solved  involve  the  fourth  power  of  r.  Wynn  [1]  demonstrated  by  numerical  methods  that, 
using  the  five  independent  magnetic  gradient  measurements  only,  there  were  4  solutions  of 
the  quadratic  equations  for  bearing  vector  and  scaled  moment  vector,  with  one  being  the 
physical  solution,  the  second  being  called  the  ‘ghost’  solution,  and  two  more  solutions 
obtained  by  reflecting  the  first  two  through  the  origin  centered  on  the  measurement  point.  It 
was  also  demonstrated  by  Wynn  [2]  that,  if  an  additional  five  measured  quantities  are 
provided  given  by  the  rate  of  change  of  the  gradient  tensor,  G'-  =  t/G(;  /  dt ,  the  solution 

represented  by  the  three  scaled  vectors  (moment,  position,  and  velocity,  W  =  V  /  r )  is  unique 
except  for  the  reflection  through  the  origin.  This  degeneracy  can  easily  be  dealt  with  in 
practical  situations. 

If  the  magnetic  field  measurement  is  usable,  we  obtain  three  additional  measurements 
for  a  total  of  eight,  allowing  us  to  resolve  both  the  range-moment  ambiguity  and  the  ghost 
problem. 
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Despite  the  difficulties,  there  are  advantages  in  using  magnetic  gradiometry  for 
locating  targets.  One  key  advantage  is  the  suppression  of  many  sources  of  noise.  The  gradient 
is  obtained  by  spatial  differentiation  of  the  magnetic  field.  Thus,  while  the  field  from  a  dipole 
decays  as  1/r3,  its  gradient  decays  as  1/r4.  This  extra  power  makes  the  effect  of  distant  noise 
sources  diminish  rapidly.  Of  course,  the  signal  from  the  target  decays  with  the  same 
functional  form,  but  gradiometry  can  provide  a  net  increase  in  signal  to  noise  ratio,  and  hence 
a  net  increase  in  detection  range. 

In  addition  to  noise  reduction,  tensor  gradiometry  provides  a  useful  quantity  that  has  a 
unique  property  in  the  sense  that  it  varies  directly  with  the  proximity  of  the  magnetic  dipole. 
Such  quantity  must  be  invariant  to  a  rotation  of  the  coordinate  system,  but  this  necessary 
condition  is  not  sufficient.  For  example,  the  scalar  representing  the  total  field  magnitude: 


B  =  B  = 


B 


x 


+  B]  +  B\ 


11/2  _  //0  111 

L, 

(m,r) 

h3  i 

m  r 

(3) 


is  invariant  to  rotation,  but  it  can  either  grow  or  shrink  with  proximity,  depending  on  the 
relative  orientation  of  the  target  dipole,  and  the  measurement  location. 

Target  localization  using  the  magnetic  field  is  the  traditional  method  used,  for 
example,  in  underwater  mine  algorithms.  The  moving  target  is  at  the  closest  point  of  approach 
(CPA)  to  the  observation  point  when  the  position  and  velocity  vectors  are  perpendicular, 
<r,V>  =  0.  It  is  known  that  the  point  of  maximum  magnetic  flux  density  is  not  directly  related 
to  the  closest  position  of  the  target.  This  can  be  easily  seen  by  taken  the  derivative  of  equation 
(3)  with  respect  to  time  and  equating  to  zero  to  obtain  the  time  at  which  the  magnetic  flux 
density  attains  its  peak  value: 


flf|B 

dt 


—  A^t  T  A^t  +  A^t  T  +  A^t  T  A ^  —  0 


(4) 


The  solution  of  this  equation  consists  of  one  real  root  and  two  conjugate  pairs  of  complex 
roots.  The  real  root  is  a  complicated  function  of  the  three  angles  between  the  vectors  r 
(arbitrary  point  on  target  trajectory),  V  and  m: 
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a  =  (2  cos2  (pn,  +  4  cos2  <prm  -  cos2  cpvm  +  6  cos  cprv  cos  (prm  cos  (pvm  + 1)/  d 
b  =  ( 2  cos (pn,  cos2  +  7 cos cos ^vm  +  3 cos (pn )/ d 

c  =  (4cos<pn  cos2  (pm  —  cos cpnn  cos (p vm  +  cos (prv )/ d 
d  =  3  cos 2  <pvm  +  1 
e  =  V2  ^3 

/  =  -a2b2  +  4b3c  +  4a3  -ISabc  +  27c2 
g  =  (9^/7  -  V3(2Z?3  -  9 ab  +  27c))'  3 
?7  =  |v|/|r| 

One  can  see  that,  by  letting  r  =  rCpA,  only  the  terms  in  cos(cprv)  equal  zero  so  that  the 
peak  time  (5)  is  not  zero  in  the  general  case,  and  the  corresponding  value  of  the  magnitude  of 
the  magnetic  flux  is  not  the  maximum  value. 

Rotation  invariant  quantities  associated  with  the  gradient  tensor  become  apparent 
when  one  attempts  to  find  its  eigenvalues  and  eigenvectors.  Using  standard  procedures, 
Pedersen  and  Rasmussen  [3]  have  shown  that  the  invariants  are  given  by: 


T  =  trace(G )  =  V  •  B  =  0 

Q  =  -Urace{ G2  )=  (|r|2  +2  rz2  )  (?) 

D  =  det(G)  =  -^^rz  (|if  +r72) 

lrl 

Both  Q  and  D  require  all  5  independent  components  of  the  tensor  G.  Q  has  only 
quadratic  as  opposed  to  cubic  terms.  The  above  equations  are  valid  in  a  coordinate  system 
obtained  after  a  rotation  such  that  a  magnetic  dipole  of  arbitrary  orientation  in  the  original 
reference  frame  becomes  aligned  with  the  Z-axis  of  the  new  reference  frame.  The  quantity  Q 
depends  directly  on  the  proximity  of  the  magnetic  object  and  it  is  often  used  to  determine 
detection  thresholds  and  to  asses  performance  of  the  system.  It  is  important  to  know  these 
quantities  because  the  measurement  platforms  are  seldom  stable. 

In  comparison  to  equations  (4-6),  the  rotation  invariant  quantity  Q  in  equation  (7)  is  a 
simple  function  of  target  position.  One  can  see  that  both  quantities,  D  and  Q,  have  the  minima 
at  |r|  =  rz,  i.e.  rx  =  rY  =  0,  but  D  is  zero  along  the  lines  rz  =  0.  Equations  (7)  were  written  in  a 
rotated  coordinate  system  where  the  new  Z-axis  coincides  with  the  direction  of  the  target 
magnetic  moment.  If  the  target  position  in  the  old  system  is  given  by  the  vector  r,  in  the  new 
system: 


4=4  +  rY+rz 


(8) 
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and,  for  a  target  moving  in  the  horizontal  plane  (rz  =  constant),  the  minimum  of  Q  occurs  at 
the  minimum  value  of  r\  +  rY  =  r(2rA ,  so  that  it  corresponds  exactly  to  the  CPA  position  of 
the  target  relative  to  the  observation  point. 

Another  quantity  often  used  in  magnetic  data  interpretation  is  the  analytic  signal, 
which  is  a  vector  encompassing  the  horizontal  derivatives  of  a  harmonic  quantity  (with  zero 
Laplacian)  and  their  Hilbert  transform.  The  amplitude  of  the  analytic  signal  can  be  used  as  an 
edge-detection  tool,  particularly  when  the  magnetic  sources  of  interest  are  relatively  close  to 
the  sensor.  Without  entering  into  the  details  [4],  one  can  define  the  total  field  magnitude  (3)  as 
the  amplitude  of  the  analytic  signal  derived  from  the  magnetic  potential  field.  In  a  similar 
way,  the  amplitude  of  the  analytic  signal  derived  from  the  magnetic  field  (B)  is  the 
Pythagorean  sum  of  the  nine  components  of  the  gradient  tensor  (tensor  magnitude): 


G  = 


P  XX  +  G  XY  +  Gjft  +  Gyx  +  Gyy  +  GyZ  +  Gyx  +  GZY  +  Gzz 


(9) 


In  summary,  there  are  specific  advantages  of  using  the  magnetic  tensor  gradiometry 
for  target  localization: 

1 .  Common  mode  rejection  of  geomagnetic  variations. 

2.  Redundancy  of  tensor  components  gives  inherent  error  correction  and  noise  estimates. 

3.  Allows  direct  determination  of  3-D  analytic  signal  (defines  source  outlines). 

4.  Rotationally  invariant  quantities  exist  that  appear  to  have  higher  resolving  power  than 
the  analytic  signal  [3]  because  of  their  faster  fall-off  rate  with  distance. 
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3.  Tracking  algorithm 


3.1  Recursive  Bayesian  estimation 


The  tracking  problem  requires  estimation  of  the  state  vector  (target  co-ordinates, 
velocity,  and  magnetic  moment)  of  a  system  that  changes  over  time  using  a  sequence  of  noisy 
measurements  (observations)  made  on  the  system.  For  the  specific  application  regarding  this 
study,  the  target  dynamics  (the  system  model)  is  described  by  a  linear  equation,  f(»),  while  the 
system  observation  (the  measurement  model)  equation,  h(»),  is  highly  non-linear.  We  assume 
that  these  models  are  available  in  a  probabilistic  form: 

P(xk  I  xa-i)  :  xa  =  f  (Vi.vt-i)  ( 1 0) 

P(z  aIxa):  z*  =h(x*,wt)  (11) 

where  xk  is  the  nx-dimensional  state  vector  of  the  system  at  time  step  k,  zk  is  the  nz- 
dimensional  observation  vector,  and  vk  and  wk  are  vectors  representing  the  process  and 
measurement  noise,  respectively.  They  have  the  dimensions  nv  and  nw.  It  is  assumed  that  the 
noise  vectors  are  independent,  identically  distributed  (i.i.d.)  random  samples  and  they  are 
independent  of  current  and  past  states. 

From  the  Bayesian  perspective,  it  is  required  to  estimate  p(xk  |  zkk)  assuming  that  the 
pdf  at  time  (k-1),  p(xk_!  |  z1:k_!),  is  available.  The  first  step  in  this  process  is  called  prediction 
and  makes  use  of  equation  (10),  which  is  assumed  to  describe  a  Markov  process  of  order  one: 

P(*k  lzl*-l)  =  jV(XA  lx*-i)/>(x*-i  \zu-i)dxk-i  (12) 


The  second  step,  the  measurement  update,  uses  the  most  recent  observation  to  produce 
the  desired  pdf  via  Bayes’  rule: 


I*™) 

P(Zk  I  Z1:A-1 ) 

P(Zk  I  Zl:A-l)  =  j P(Zk  I  XA-)  P(Xk  I  Zl:A-l)  dxk 


(13) 


where  the  second  equation  is  the  normalization  constant.  Once  the  posterior  pdf  is  determined, 
it  is  straightforward  conceptually  to  produce  any  desired  statistic  of  xk.  For  instance,  the 
minimum  mean-square  error  (MMSE)  estimate  of  the  current  state  could  be  found  by 
computing  the  conditional  mean: 


x 


k 


jXA  P(*k  I  Z\:k  )dxk 


(14) 


The  conditional  covariance  matrix  is  obtained  in  a  similar  way: 
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P1  =  J(XA  -XA-f  P(*k\Z\±)d*k 


(15) 


In  general,  the  recursive  propagation  of  the  posterior  density  cannot  be  determined 
analytically  because  the  integrals  in  (12)  and  (13)  do  not  have  closed-form  solutions. 
Solutions  do  exist  in  a  restrictive  set  of  cases.  For  example,  if  f(»)  and  h(»)  are  linear  functions 
and  if  Gaussian  distributions  are  assumed  for  x,  v,  and  w,  the  estimation  of  states  is  reduced  to 
the  well-known  Kalman  filter. 

The  problem  of  tracking  a  magnetic  dipole  does  not  satisfy  the  original  Kalman  filter 
requirements  because  the  system  observation  is  non-linear.  Moreover,  because  the  target  can 
approach  the  sensors  from  any  direction  and  can  maneuver  at  any  time,  the  true  posterior 
density  is  multi-modal  and  a  Gaussian  description  will  be  inaccurate. 

In  order  to  deal  with  non-linear  systems  and/or  non-Gaussian  reality,  two  categories  of 
techniques  have  been  developed:  parametric  and  non-parametric.  The  parametric  techniques 
are  based  on  improvements  of  the  Kalman  filter.  These  filters  (for  example,  extended  and 
unscented  Kalman  [5]  filters)  can  handle  non-linear  equations,  but  they  implicitly 
approximate  the  posterior  density  as  Gaussian.  The  non-parametric  techniques  are  based  on 
Monte  Carlo  simulations  and  are  the  subject  of  the  present  study.  These  filters  assume  no 
functional  form,  but  instead  use  a  set  of  random  samples  (particles)  to  estimate  the  posteriors. 
The  advantage  is  that  the  particle  filters  can  accommodate  simultaneous  alternative 
hypotheses  that  can  describe  a  multi-modal  distribution  well. 


3.2  Particle  filter  implementation 


The  basic  idea  of  the  Monte  Carlo  based  approach  to  an  intractable  Bayesian  filtering 
case  is  to  approximate  an  unknown  distribution,  p,  by  a  set  of  properly  weighted  particles 
drawn  from  a  known  distribution  q.  In  this  way,  the  difficult  problem  of  distribution 
estimation  is  converted  to  an  easy  problem  of  weight  estimation.  The  exact  form  of  the 
proposal  distribution  q  is  a  critical  issue  in  designing  the  particle  filter  and  is  usually 
approximated  to  facilitate  easy  sampling. 

A  numerical  approximation  to  the  recursive  Bayesian  filtering  method  given  by  the 
equations  (12)  and  (13)  is  the  following  algorithm  [6]: 


1.  Initialization:  sample  N  particles  xk(l),  i  =  1,  2,  ...,  N,  from  the  proposal  distribution. 
The  proposal  distribution  can  be  the  transition  prior  as  used  in  the  conventional  particle 
filters,  or  more  advanced  distributions  like  the  one  used  in  this  study. 

2.  Measurement  update:  update  the  importance  weights.  The  Bayesian  sequential 
importance  sampling  (SIS)  procedure  gives  a  recursive  calculation  of  the  normalized 
weight  [6]: 
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As  an  approximation  to  (14)  take: 


i=l,N 


X 


(0 

k 


3.  Re-sampling  is  a  necessary  step  introduced  in  particle  filtering  algorithms  to  reduce 
the  degeneration  of  samples.  In  practice  it  was  noticed  that,  after  a  few  iterations,  one  of 
the  importance  weights  tends  to  one,  while  the  others  become  zero.  To  avoid  the 
degeneracy,  the  sampling  importance  re-sampling  (SIR)  method  selects  N  samples  with 
replacement  from  the  set  xkw,  where  the  probability  to  take  sample  ‘i’  is  wk(l>.  Then  set 
wk(i)=  1/N,  i=  1,2,  ...,N. 

4.  Prediction:  assuming  that  the  probability  of  the  process  noise  is  known,  use  equation 
(6)  to  simulate  xk+1(l),  i  =  1,2,  . . .,  N. 

5.  Set  k  =  k  +1,  and  iterate  to  item  2. 


3.3  The  unscented  Kalman  filter 


As  mentioned,  the  deficiency  of  the  sequential  importance  sampling  (SIS) 
approximation  is  that  the  proposal  distribution  may  be  very  different  from  the  posterior 
distribution,  especially  if  using  the  transition  prior  as  the  proposal  distribution.  An  improved 
proposal  distribution  must  incorporate  the  current  observation  data  with  the  optimal  Gaussian 
approximation  of  the  state. 

In  a  previous  study  [7]  on  the  magnetic  dipole  tracking  application,  it  was  shown  that 
the  unscented  Kalman  filter  (UKF)  is  the  best  Kalman  filter  for  the  non-linear  systems.  The 
UKF  is  so  named  because  it  implements  the  Kalman  recursion  using  the  sample  points 
provided  by  the  unscented  transform.  The  unscented  transform  deterministically  generates  a 
set  of  points  that  have  a  certain  mean  and  sample  covariance.  The  non-linear  function  is  then 
applied  to  each  of  the  sample  points,  yielding  a  transformed  sample  from  which  the  predicted 
mean  and  covariance  are  calculated.  The  estimate  of  the  conditional  mean  provided  by  the 
UKF  is  shown  to  be  correct  up  to  the  second  order  of  its  Taylor  series  expansion.  Reference 
[5]  gives  the  implementation  of  UKF  algorithm,  which  is  presented  in  the  Annex. 

Because  the  UKF  is  the  best  to  accurately  propagate  the  mean  and  covariance  of  the 
Gaussian  approximation  to  the  state  distribution,  it  can  be  used  to  generate  the  proposal 
distribution  for  the  particle  filter.  In  this  way,  one  obtains  a  parametric/non-parametric  hybrid 
filter  called  the  unscented  particle  filter  (UPF). 
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4.  Modeling  the  target  dynamics 


For  a  full  characterization  of  the  target,  the  entire  system  at  time  step  k  can  be 
represented  by  the  state  vector: 

xa-  =  irx  rY  r?  Vx  VY  mx  mY  mz)T  (17) 

The  only  simplifying  assumption  made  in  equation  (17)  is  that  the  target  is  moving 
horizontally. 

As  mentioned  before,  the  recursive  Bayesian  estimation  method  is  based  on  the 
assumption  that  the  state  variables  follow  a  first  order  Markov  process.  Thus,  the  discrete 
equations  of  the  state  evolution  of  the  target  are  obtained  using  the  piece-wise  approximation: 

rx  (*)  =  rx (k-l)  +  At  Vx (k  - 1) 

Vx{k)  =  Vx{k- 1)  (18) 

mx(k)  =  mx[k- 1) 


where  At  is  the  time  increment  between  the  data  samples  in  seconds,  and  similar  relations 
exist  for  the  Y  and  Z  components  (Vz  =  0). 

In  addition  to  the  robustness  against  the  noise  in  data,  the  proposed  statistical 
recursive  method  has  another  advantage  in  comparison  to  the  direct  inversion  in  [1].  We  know 
[1]  that,  using  the  five  independent  magnetic  gradient  measurements  only,  there  is  a  ‘ghost’ 
solution  for  bearing  vector  and  scaled  moment  vector  in  addition  to  the  physical  solution.  Two 
more  solutions  are  obtained  by  reflecting  the  first  two  through  the  origin  centered  on  the 
measurement  point.  In  this  paper,  we  demonstrate  that  the  imposed  evolution  of  the  state 
vector  in  equation  (18)  can  resolve  the  ‘ghost’  ambiguity  of  solutions  while  using  only 
magnetic  gradient  information.  Let  us  assume  that  at  time  step  (k- 1 )  the  state  vector  does  not 
represent  the  physical  scaled  solution,  R,  W  and  M,  but  the  ‘ghost’  solution.  Even  so,  the 
estimated  value  of  the  gradient  at  time  (k-1)  would  correspond  to  the  measured  value.  Using 
equations  (18)  and  the  Taylor’s  series  expansion,  one  can  write  to  the  first  approximation  the 
gradient  at  time  k: 


Gy  (R(A-) ) -  Gy  (R(a-d  +  W(A_  1)  A?) ~  G,  (R(a_i)  )  + 


5G„. 


5R 


W(*_n  At  = 


dt 


(19) 


«=*-!) 


At 


As  shown  before  [2],  the  ‘ghost’  solution  can  produce  the  correct  gradient,  but  it  cannot 
produce  the  correct  gradient  rate,  irrespective  of  the  value  of  the  scaled  velocity. 
Consequently,  from  equation  (19),  the  gradient  value  at  time  (t  =  k)  will  differ  from  the 
measured  value.  This  means  that  the  proposed  state  transitions  and  measurement  model  (five 
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gradient  components),  to  be  able  to  reproduce  the  measurements,  must  accept  only  the 
physical  solution  (and  the  reflection).  It  is  easily  seen  in  equation  (2)  that,  if  m  and  r  are 
replaced  by  -m  and  -r,  respectively,  the  value  of  G  does  not  change.  Necessarily,  the  scaled 
velocity  WfyM)  will  produce  the  observed  gradient  rates. 

To  obtain  an  un-ambiguous  magnetic  tracking  one  needs  three  additional  independent 
measurements  of  the  magnetic  field  for  a  total  of  eight  data  points  [8].  The  use  of  a  stationary 
gradiometer  and  magnetometer  resolves  both  the  ‘ghost’  problem  and  the  range-moment 
ambiguity.  In  this  case,  the  solution  is  unique  because  the  replacement  of  m  and  r  by  -m  and 
-r,  respectively,  in  equation  (1)  changes  the  sign  of  B. 

In  the  following,  both  these  two  possibilities  will  be  analyzed.  In  the  first  case,  the 
measurement  vector  at  time  k  has  only  gradiometer  data: 

zk  =  (G xx  GYX  G xz  Gxy  Gyy  Gzy  Cxz  Gyz  Gzz  )  (20. a) 


Only  five  components  are  independent,  but  because  the  data  acquisition  chanels  may 
not  be  identical  all  nine  components  will  be  used  in  the  localization.  Gradients  alone  will  give 
the  bearing  vector  to  the  target  and  the  orientation  of  its  moment. 

For  the  second  case,  the  measurement  vector  will  be  a  combination  of  magnetic 
gradient  and  magnetic  field  components: 


Zk  ~  (G  XX  Gyx  GXZ  Gyy  Gyy  Gyy  By,  By  By,  ) 


(20.b) 


This  approach  resolves  the  range-moment  ambiguity  and  eliminates  ‘ghosts’. 

The  equations  (17)  and  (18)  are  the  process  and  measurement  equations,  respectively. 
As  one  can  see,  the  process  function  !'(•)  in  equation  (14)  is  linear,  and  the  measurement 
function  h(»)  in  equations  (20),  (2)  and  (1)  is  highly  non-linear. 
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5.  Experimental  results 

Having  outlined  the  background  to  the  problem  and  suggested  how  the  tracking  of  a 
magnetic  dipole  target  in  real  time  might  be  carried  out,  these  ideas  are  next  applied  to  a  set  of 
real  data  collected  with  a  gradiometer  sensor.  The  experiment  took  place  in  Kejimkujik 
National  Park,  Nova  Scotia,  where  the  ambient  magnetic  activity  due  to  external  sources,  such 
as  power  lines  and  traffic,  is  very  low.  On  the  satellite  image  of  the  experiment  area  (Fig.l), 
the  cross  indicates  the  position  of  the  gradiometer  and  the  double-headed  parallel  arrows 
approximately  represent  the  tracks  used  in  the  experiment.  Each  track  has  a  number  attached 
that  represents  the  perpendicular  distance  (CPA)  in  meters  from  the  gradiometer  to  that  track. 
The  target  was  an  automobile  that  was  driven  on  the  marked  (5,  10,  14,  40  and  45  meters) 
staight  paths  towards  North-West  and  South-East,  and  on  several  curved  paths  not  shown  in 
the  picture. 

The  gradiometer  sensor  is  placed  at  the  center  of  a  Cartesian  coordinate  system  with 
the  X-axis  oriented  towards  the  magnetic  North,  Y-axis  towards  West,  and  the  Z-axis  oriented 
upwards  with  zero  at  the  ground  level.  Consequently,  the  target  coordinates  will  be  given  in 
this  reference  frame.  Note  that  on  the  satellite  picture  the  tracks  are  drawn  in  the  geographical 
coordinates.  A  yellow  arrow  in  Fig.l  with  the  label  ‘Nm’  approximately  indicates  the 
magnetic  North  in  the  experiment  area  (where  the  declination  is  21°),  so  that  the  actual  angle 
between  the  X-axis  and  the  staight  tracks  is  about  45°. 

The  gradiometer  sensor  was  built  with  4  tri-axial  Bartington  fluxgate  magnetometers, 
one  placed  in  the  origin,  and  the  other  three  on  the  orthogonal  axes  at  0.65m  from  origin. 
Estimates  of  the  gradient  tensor  components  are  made  by  forward  differencing  the 
corresponding  magnetic  field  components.  The  three  independent  magnetic  field  components 
used  in  this  study  additional  to  the  gradiometer  data  are  measured  at  the  magnetometer  sensor 
placed  in  the  origin.  All  data  were  acquired  at  a  sampling  rate  of  10Hz.  The  maximum 
measured  noise  level  of  this  system  was  detected  in  the  GZz  channel  and  was  90  pT/m 
rms/Vl  Iz  in  the  0-5Hz  bandpass.  The  individual  sensor  noise  level  measured  at  the  central 
magnetometer  was  40  pT  rms/Vl  Iz  (the  specification  for  the  Bartington  magnetometer  alone  is 
<12  pT  rms/N  I  Iz  at  1Hz).  The  measured  noise  depends  on  both  the  gradiometer  construction 
[8]  and  the  instrumental  noise.  Due  to  this  relatively  high  level  of  noise,  this  system  exhibits 
rather  restricted  range,  but  it  provided  a  full  test  for  the  tracking  algorithm  against  a  real 
dipole  source.  For  comparison,  the  superconducting  magnetic  gradiometer  used  for  the 
tracking  experiment  in  [8]  had  a  noise  level  of  about  0.1  pT/m  rms/Vl  Iz  at  1Hz. 

In  applying  the  Bayesian  filtering  technique  to  the  system,  the  initial  conditions  and 
the  noise  covariance  matrixes  need  to  be  specified.  In  the  initialization  step,  the  particles 
should  be  drawn  from  an  unknown  proposal  distribution.  The  basic  assumption  is  that  the 
target  can  approach  the  sensors  from  any  horizontal  direction.  Therefore,  the  filter  must 
accommodate  simultaneous  alternative  hypotheses  until  they  can  remove  the  ambiguity  by 
future  measurements.  A  reasonable  initial  estimate  of  the  horizontal  position  is  an 
approximate  circle  around  the  gradiometer  sensor  with  a  radius  of  about  100m  from  where  the 
magnetic  signal  becomes  sizable.  In  the  present  study,  36  particles  were  used  with  the 
horizontal  positions  spread  over  a  circle  every  10°  from  0°  to  350°.  Because  we  have  no  a 
priori  information  about  the  vertical  position,  the  magnitude  of  velocity,  and  magnetic  dipole 
moments,  a  good  initial  estimate  of  these  parameters  is  merely  the  null  vector. 
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The  initial  covariance  matrix,  P(0|0),  gives  a  measure  of  belief  in  the  initial  state 
estimate.  It  is  assumed  that  initially  all  the  states  are  un-correlated,  so  that  the  matrix  is 
diagonal.  This  matrix  is  not  known  and  has  to  be  sufficiently  large,  but  the  initial  P(0|0)  is 
forgotten  as  more  data  is  processed.  The  measurement  noise  covariance  matrix  can  be 
estimated  directly  from  the  actual  data  and,  once  calculated,  it  does  not  change  during  the 
filter  run. 

The  process  noise  covariance  is  zero  for  a  deterministic  process.  However,  it  was 
practically  proved  to  be  a  good  idea  to  introduce  random  perturbations  in  the  target  position 
and  velocity.  These  small  perturbations  account  for  the  target  maneuvers  and  prevent 
divergence,  so  that  the  process  noise  covariance  may  be  regarded  as  a  tuning  parameter  of  the 
filter. 

To  demonstrate  the  tracking  algorithm,  we  use  the  data  collected  with  the  target 
moving  straight  towards  North-West  on  the  path  situated  at  10  meters  CPA  from  the  sensor. 
All  9  components  of  the  magnetic  tensor  gradient  are  plotted  in  Fig.2  in  nano-T/m.  Due  to  the 
instrumentation  errors  the  symmetry  of  the  tensor  is  not  perfectly  preserved  (see,  for  example, 
that  Gxz  is  slightly  different  from  GZx)-  The  geometry  of  the  experiment  is  shown  in  Fig.  3  in 
the  sensor  coordinate  system  with  X-axis  towards  the  magnetic  North.  The  first  case  analyzed 
is  when  the  measurement  vector  has  only  gradiometer  data  (20. a).  Figure  4  shows  the  bearing 
vector  that  proves  to  follow  very  well  (dotted  line)  the  target  in  motion.  As  discussed  above, 
this  solution  is  not  unique  and  the  filtering  method  could  produce  as  well  the  reflection  of  the 
bearing  vector  through  the  origin.  The  other  state  variables,  velocity  and  magnetic  moment, 
are  obtained  as  scaled  parameters,  so  that  they  cannot  characterize  the  target. 

A  unique  solution  in  absolute  units  is  obtained  when,  in  addition  to  the  5  independent 
gradient  components,  the  magnetic  field  components  measured  at  the  magnetometer  placed  in 
the  origin  of  the  coordinate  system  are  included  in  the  observation  vector  (20. b).  Figures  5  to 
8  present  the  evolution  in  time  of  several  system  state  variables  that  characterize  the  target  in 
motion.  One  can  see  that  the  trajectory  of  the  target  follows  very  closely  the  10m  track 
represented  by  the  blue  arrow  in  Fig.5.  The  vertical  position  of  target  is  between  zero  and  5m. 
Other  parameters,  such  as  the  target  velocity  (~  3m/sec)  and  Z-magnetic  moment  (-0.5kA-m2, 
oriented  downwards  as  expected  at  this  latitude),  are  not  known  a  priori  in  this  experiment, 
but  their  estimated  values  are  reasonable  for  this  type  of  target.  The  other  two  components  of 
the  magnetic  moment  are  practically  zero. 

Finally,  an  interesting  test  of  the  algorithm  makes  use  of  data  obtained  for  a  curved 
track.  The  test  shows  the  capability  of  the  algorithm  to  adapt  to  the  target  maneuvers.  Even  if 
the  target  motion  is  assumed  rectilinear  in  equations  ( 1 8),  the  small  perturbations  introduced 
by  the  process  noise  allows  for  a  variable  velocity  which,  in  turn,  modify  the  trajectory  from 
rectilinear  to  curve.  If  the  trajectory  is  curved,  one  expects  modifications  in  the  behavior  of 
the  rotationally  invariant  quantities  Q  and  B.  For  a  rectilinear  movement  there  is  one  CPA 
point  so  that  Q  has  only  one  minima  (-Q  and  B  have  a  maxima)  presented  in  Fig.9  using  the 
previous  data  recorded  with  the  target  moving  straight  on  the  10m  track.  Also  Fig.9  presents 
the  estimated  inverse  value  of  the  target  position  vector  magnitude  for  comparison.  There  is  a 
delay  of  about  0.2  sec  (2  data  samples)  between  the  CPA  time  determined  by  Q  (or  B)  and  the 
one  estimated  by  1/r,  possibly  because  the  filter  requires  a  certain  amount  of  data  to  correct 
the  system  dynamics. 

Anticipating  a  bit  the  results,  the  same  quantities,  -Q,  B  and  estimates  of  1/r,  are 
plotted  in  Fig.  10  for  the  target  moving  on  the  curved  path.  The  target  came  closer  to  sensor 
twice,  but  not  at  the  same  distance.  The  CPA  point  given  by  the  estimated  values  of  1/r 
corresponds  to  the  highest  peak  of-Q.  This  is  not  the  case  for  the  maximum  of  B.  The  plot 
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clearly  indicates  that  the  quantity  Q  accurately  determines  the  variations  in  the  target 
proximity. 

The  rest  of  figures  from  11  to  14  present  the  most  important  parameters  of  the  target 
when  moving  on  the  curved  track  as  estimated  by  the  unscented  particle  filter.  The  blue  arrow 
in  Fig.  11  indicates  the  straight  10m  path.  The  target  trajectory  in  Fig.  11  is  assumed  to  be 
correctly  estimated  because  there  is  no  tracking  information  collected  during  this  experiment. 
Flowever,  the  target  was  detected  to  move  on  the  z  =  0  horizontal  plane  (Fig.  12)  and  can  be 
characterized  by  its  Z-component  of  magnetic  moment  that  was  estimated,  as  previously,  to 
be  -0.5kA-m2. 

Similar  results  as  above  were  obtained  for  the  target  moving  on  the  5  and  14  meters 
tracks.  The  value  of  the  magnetic  moment  explains  why  there  was  no  signal  recorded  by  the 
gradiometer  when  the  automobile  moved  on  the  40  and  45  meters  tracks.  Considering  the 
decay  with  the  fourth  power  of  the  distance,  the  magnitude  of  the  signal  when  the  target  is  at 
40m  from  the  sensor  would  be  of  about  20  pT/m,  which  is  well  below  the  level  of  the 
environmental  noise  for  this  system. 
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Figure  1.  Satellite  view  of  the  Keji  park  layout  of  tracking  experiment. 
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Figure  2.  All  9  magnetic  gradient  components  for  the  target  on  the  1 0m  track. 
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Figure  3.  Experiment  coordinate  system  and  the  1 0m  track. 


Figure  4.  The  estimated  and  true  bearing  vector  of  target  on  the  10m  track. 
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Target  Z  position  (m)  X-axis  (to  North) 


Figure  5.  X-Y  estimates  of  target  position  (10m  path). 


Figure  6.  Z  estimates  of  target  position  (10m  path). 
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Figure  7.  Estimated  target  Z-magnetic  moment  component  (10m  path). 


Figure  8.  Estimated  magnitude  of  target  velocity  (10m  path). 


DRDC  Atlantic  TM  2006-188 


Figure  9.  Plot  of-Q,  B  and  estimated  103/\r\  for  the  target  on  10m  track. 


Figure  10.  Plot  of-Q ,  B  and  estimated  103/\r\  for  the  target  on  the  curved  track. 
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6.  Conclusions 


This  report  investigates  the  possibility  of  using  the  unscented  particle  filter  (UPF)  for 
tracking  and  classification  of  a  magnetic  dipole  with  a  gradiometer.  A  fluxgate  gradiometer- 
magnetometer  array  was  constructed  and  used  to  demonstrate  the  feasibility  of  the  processing 
technique  against  a  real  target.  It  was  shown  that  the  effective  range  of  operation  for  this 
system  is  given  by  the  actual  noise  level  that  depends  on  gradiometer  construction  and 
instrumental  noise. 

The  magnetic  gradient  measurements  offer  some  advantages  over  the  magnetic  field. 
In  addition  to  noise  reduction,  tensor  gradiometry  provides  a  useful  quantity  (Q)  that  has  a 
unique  property  in  the  sense  that  it  varies  directly  with  the  proximity  of  the  magnetic  dipole. 
This  quantity  is  invariant  to  a  rotation  of  the  coordinate  system,  which  is  another  important 
property  because  the  measurement  platforms  are  seldom  stable.  This  report  shows  that  the 
quantity  Q  appears  to  have  higher  resolving  power  than  the  magnitude  of  the  magnetic  field, 
which  is  also  rotationally  invariant. 

The  report  investigates  the  inverse  problem  of  the  determination  of  the  position, 
velocity  and  magnetic  parameters  of  the  target  at  a  given  time  step  from  its  magnetic  signature 
collected  up  to  and  including  that  time  step.  Two  cases  are  separately  analyzed:  (i)  the 
observation  contains  only  gradiometer  data  when  a  double  solution  given  by  scaled  state 
vectors  exists,  and  (ii)  magnetic  field  components  are  added  to  the  previous  case  and  all  target 
parameters  are  uniquely  determined. 

In  the  first  case,  the  determination  of  the  azimuth  and  the  elevation  of  the  target 
relative  to  the  gradiometer  sensor  represents  substantially  more  information  than  provided  by 
a  vector  magnetometer.  However,  the  solution  is  a  scaled  vector  (range-velocity-moment 
ambiguity)  and  it  is  accompanied  by  its  reflection  through  the  origin. 

The  situation  is  greatly  improved  if  the  magnetic  field  measurement  is  usable.  In  this 
case,  one  obtains  three  additional  independent  measurements  for  a  total  of  eight,  resolving 
both  the  ambiguity  and  the  solution  uniqueness  problem.  Compared  to  the  non-recursive 
solution  of  the  dipole  field  equations  yielding  position  and  dipole  moment  due  to  Wynn  [2,  8], 
the  present  statistical  recursive  method  is  easier  to  implement  and  is  robust  against  the  noise 
in  the  data,  which  may  be  a  serious  problem  when  direct  inversion  is  used. 

The  filtering  technique  applied  on  the  experimental  data  allows  all  the  target 
parameters  to  be  estimated  with  a  good  level  of  accuracy.  The  proposed  method  of  recursive 
Bayesian  estimation  proved  to  be  applicable  to  various  situations  (straight  or  curved  tracks). 
The  validation  of  the  presented  algorithms  on  experimental  data  indicates  the  possibility  to 
incorporate  it  into  automated  surveillance  systems. 
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Figure  11.  X-Y  estimates  of  target  position  (curved  path). 
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Figure  12.  Z  estimates  of  target  position  (curved  path). 
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Figure  13.  Estimated  target  Z-magnetic  moment  component  (curved  path). 
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Figure  14.  Estimated  magnitude  of  target  velocity  (curved  path). 
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Annexes 


Implementing  the  unscented  Kalman  filter 

The  unscented  transformation  uses  a  set  of  weighted  samples  or  ‘sigma  points’,  Sj  =  {W,,  %i}, 
to  completely  capture  the  true  mean  and  covariance  of  the  random  variable  x  and  then 
propagates  the  sigma  points  through  the  non-linear  function.  We  use  the  following  notation:  K 
is  the  Kalman  gain,  W;  are  the  weights,  /.  is  the  scaling  parameter,  na=nx+nv+nw,  and 

xfl  =[xr  vr  wrf,  Ka  =  (xx)r  (Kv)r  (Kw)r  r 

1 .  Initialize  with: 

x0=£[x0];  P0=^[(x0-x0)(x0-x0fJ  x;  =£[xa]=[xo  Oof 

rp0oo_ 

P0a=i<[(xa-x0a)(x;-x0a)r]=  0Q0 

[oOR 


2.  For  discrete  time  samples,  k  =  1,2,  . . K, 
(a)  Calculate  sigma  points: 


(b)  Time  update: 

^>k|k-l  =  f^k-l^k-l) 

2 na 

X  _  V1  yy(m)  -v 

Ak|k-1  ~  Z-j  ”i 

i= 0 

f*k|k-l  =  '  [^>i,k|k-l  “  Xk|k  1  ][^i,k|k-l  _  Xk|k-1 

i=0 
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2 nc, 

*k,k-i=E^(w) 


(c)  Measurement  update  equations: 

=  *  \^i,k\k-\  ~  Zk|k  -1  ]  [Zi,t|it-I  ~  Zk|k-l] 

i=0 

^k\k  ~  1  —  Xk|k-l][ZUM  _  Zk|k  1  ] 


i=0 


Where  Pzz,  Pxz  are  respectively  the  covariance  matrix  of  the  measurement  and  the  cross¬ 
covariance  of  the  measurement  and  the  state  variable.  Next: 


is  _  p  XZ  pZZ 
~  rk\k  [rk\k 


Xk  _  Xk|k-1  +  Kk  (zk  Zk|k-1  ) 

P  — P  i.  p22ivT 

rk  _  rk|k-l  ^k  ri|i-  ^k 
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